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ABSTRACT 


It  is  demonstrated  how  the  kernel  functions  of  single  integral 
constitutive  equations  may  he  determined  from  analysis  of  experiments  in  time 
dependent  shear-free  flows  alone.  It  is  assumed  that  the  kernel  functions  may 
be  factored  into  a  product  of  a  linear  viscoelastic  function  and  a  finite- 
strain  dependent  function.  No  assumption  is  needed  on  the  strain-dependent 
function  except  that  it  must  be  continuous  within  the  attainable  invariant 
space.  The  experimental  data  for  ever-increasing  deformations  are  not 
compatible  with  the  assumptions  inherent  in  single  integral  constitutive 
equations  with  factorized  kernel  functions. 
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The  characterization  of  the  flow  properties  of  polymeric  ir.elts  represents 
an  important  and  still  partially  unsolved  problem.  The  "constitutive 
equations"  used  in  the  characterization  involve  unknown  kernel  functions  of 
several  variables  that  need  to  be  determined  from  experiments.  Usually  some 
specific  analytical  form  of  the  kernel  functions  with  a  small  number  of 
parameters  is  either  assumed  or  derived  from  molecular  arguments.  In  this 
report  we  show  how  one  may  determine  the  kernel  functions  of  some  rather 
general  constitutive  equations  directly  from  experiments  on  a  polymer  melt. 

The  method  is  analagous  to  the  determination  of  the  strain  energy  function  in 
crosslinked  rubber.  The  paper  should  be  useful  to  the  theoretical  rheologist 
who  is  continually  searching  for  "the  most  general"  constitutive  equation 
needed  to  characterize  polymer  melts,  and  to  the  experimental  rheologist  who 
wishes  to  analyze  measurements  with  a  minimum  of  assumptions. 
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SINGLE  INTEGRAL  CONSTITUTIVE  EQUATIONS  FOR  VISCOELASTIC  FLUIDS 


Poul  Bach  and  Ole  Hassager* 
1.  INTRODUCTION  AND  DEFINITION  OF  MODELS 


The  purpose  of  this  report  is  to  demonstrate  how  the  kernel  functions  of  two  rather 
general  single  integral  constitutive  equations  for  viscoelastic  fluids  may  be  determined 
from  experimental  measurements.  The  two  particular  equations  are: 

"The  factorized  K-BKZ  model": 


rrt,  =/  MU  -f>  l||-  If0I  ♦  Jdt. 

-CO  12 


and  "The  factorized  Rivlin  Sawyers  model": 


r(t)  =  /  M(t  -  t'lH^o,  +  ♦2It°1]df 


Here  T(t)  is  the  stress  tensor  at  a  given  particle  at  the  current  time  t  and 
[0] 

and  1  are  strain  tensors  that  differ  from  the  Finger  strain  tensor  %  and  its 
inverse  B  \  respectively  only  by  a  unit  tensor  6: 

xto]  -  i  -  S<  i[01  “  l'  ~  i 

The  components  of  g  are  given  by  Bj^  =  Ox^/Sx^)  (5x./3  x^)  where  the  x^  are  the 
coordinates  at  time  t  of  the  particle  with  coordinates  x|  at  time  t1.  In  addition  in 
Eqs.  (1.1)  and  (1.2)  the  function  M  is  the  memory  function  of  linear  viscoelasticity 
related  to  the  linear  viscoelastic  relaxation  modulus  G  as  follows 
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Also  ■<>  1  ( I  j ,  I^) ,  an<^  W(I^,I2)  are  nor>dimensional  scalar  functions  that  depend 

on  the  scalar  invariants  1^  and  I2  defined  hy 

-1 

It  =  tr  B;  I2  =  tr  B  (1.5) 

Certainly  the  model  in  Fq.  (1.1)  is  contained  in  Fa.  (1.2),  however  for  the  purpose  of  this 
report  it  is  convenient  to  consider  it  as  a  separate  model. 

In  connection  with  the  models  in  Eqs .  (1.1)  and  (1.21)  we  note  that: 

i)  If  fr'(I1,I2)  and  the  A  (x1fi  ),  i  =  1,2  are  analytic  at  (I^Ij)  =  (3,3)  then 
the  models  may  be  approximated  in  slow  flow  by  a  "third  order  fluid"  with  parameters  given 
in  Bird,  Armstrong  and  Hassager  (1985). 

ii)  Renardy  (1984)  has  derived  sufficient  conditions  on  W  under  which  initial  value 
problems  based  on  Eq.  (1.1)  are  always  well  posed,  i.e.  Hadamard  type  instability  cannot 
occur. 

iii)  Hassager  (1981)  has  derived  a  variational  principle  for  Eq.  (1.1)  useful  for 
simulation  of  creeping  flow  situations  based  on  that  constitutive  equation. 

The  model  in  Fq.  (1.1)  is  a  factorized  form  of  a  more  general  model  proposed 
independently  by  Kaye  (1962)  and  Bernstein,  Kearsley  and  Zapas  (1963): 


r,v  y 
^7*t0] 
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(1.6) 


where  V  is  a  scalar  function  of  (t  -  t'),  1^  and  I2.  We  call  this  the  K-BKZ  model. 

The  arguments  leading  to  this  model  may  be  traced  back  to  Lodge  (1964)  who  noted  the 
striking  similarity  between  the  behavior  of  polymer  melts  and  crosslinked  natural  rubber. 
Lodae  then  started  with  a  stress-strain  relation  for  crosslinked  natural  rubber  obtained  on 
the  basis  of  a  molecular  theory,  and  then  used  a  heuristic  argument  to  "liquify"  the 
relation  to  arrive  at  the  "rubberlike  liquid  model"  in  Table  1.1.  In  the  same  way  does  Eq. 
(1.6)  represent  a  heuristic  "liquif ication"  of  the  nonlinear  stress-strain  relation  for 
incompressible  isotropic  materials.  For  example  if  in  Eq.  (1.1)  one  replaces  M(t  -  t ' ) 
by  Gp*(t'  -  t^)  one  obtains  the  stress-strain  relation  for  an  elastic  material  with 


stress  free  state  equal  to  the  configuration  of  time  tp  and  with  strain  energy 


function  (GqW).  Certainly  the  special  status  given  to  one  particular  configuration  is 
incompatible  with  the  concept  of  a  fluid,  and  it  is  understood  that  the  function  V  in  Eq. 
(1.6)  or  M  in  Eq.  (1.1)  do  not  give  special  status  to  any  one  past  configuration  (but 
they  decay  to  zero  as  (t  -  t  * )  tends  to  infinity  sufficiently  fast  that  the  integrals 
converge). 

The  argumentation  leading  to  Eq.  (1.6)  was  questioned  by  Rivlin  and  Sawyers  (1971)  who 
instead  started  with  the  physical  assumption  that  the  effects  on  the  stress  at  time  t  of 
the  deformations  at  different  past  times  t'  are  independent  of  each  other,  with  this 
assumption  it  may  be  shown  from  the  theory  of  additive  functionals  (Martin  and  Mizel 
(1964),  Chacon  and  Friedman  (1965))  that  the  most  general  constitutive  equation  for 
isotropic  fluids  is: 


where  the  tit’s  are  functions  of  (t  -  t'),  1^  and  I2 .  We  call  this  the  Rivlin  Sawyers 

model . 

We  now  see  that  the  two  models  to  be  investigated  in  this  report,  Eqs.  (1.1)  and  (1.2) 
follow  from  the  assumptions  inherent  in  K-BKZ  model  and  the  Rivlin  Sawyer  model  plus  the 
additional  assumptions  that  respectively 

V(t  -  t ' , I , , I 2 )  =  M(t  -  t’)W(I1,l2)  (1.8) 

tbi  ( t  -t',1,,12)  =  M(t  -  t’  ><^<1,  ,I2>»  i  =  1/2  (1.9) 

It  is  at  this  point  reasonable  to  ask  why  the  seemingly  arbitrary  assumptions  in  Eqs.  (1.8) 
and  (1.9)  are  added  to  the  assumptions  already  present  in  the  X-BKZ  and  Rivlin  Sawyers 
models.  We  do  this  for  the  following  reasons: 

i)  Two  key  types  of  molecular  models,  the  Lodge  (1956)  network  model  on  the  one  hand 

and  the  Curtiss  Bird  (R  =  0)  (1981)  and  Doi  Edwards  (1978)  melt  models  on  the  other  hand 

show  such  a  factorization  and  may  in  fact  be  represented  by  a  K-BKZ  model  with  a  potential. 

ii)  Carefully  performed  experiments  (e.g.  Laun  (1978),  Einaga,  Osaki  and  Kurata 
(1971))  indicate  that  at  least  some  polymers  may  allow  for  such  a  factorization  in  simple 
shear  deformations. 
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iii)  From  a  practical  point  of  view  it  may  be  impossible  to  determine  a  function  of 
three  independent  variables  on  the  basis  of  a  limited  number  of  experimental  data.  For 
this  reason  a  number  of  empirical  models  with  a  limited  number  of  parameters  have  been 
proposed,  as  shown  in  Table  1.2.  These  empirical  relations  all  involve  the  assumption  in 
Ea.  (1.8). 


Table  1.1  Examples  of  factorized  K-BKZ  models 


Reference 

M(s) 

M(I,,I2) 

Adjustable  parameters 

»  M.  -s/X. 

V  i  i 

1  ~2  6 

Lodge  (19S6) 

r1 

y^.X^  for  i  ”  1,2, . . . 

i=1  X^ 

Curtiss  and  Bird3 

f-  NnkTX-1 

I 

4tt  1 

In  (B  :  uu)du 

(NnkTX  %  and  X 

(1981)  with  their 

e  =  0. 

Y  exp( -x^a2s/X ) 

a, odd 

1 

1 

Notes:  aIn  this  model  u  is  a  unit  vector,  and  j  du  is  an  integral  over  a  unit 
sphere.  This  model  is  identical  in  form  to  that  of  Doi  and  Edwards  (1978),  (1979) 


Table  1.2 


Examples  factorized  Rivlin  Sawyers  models  (empirical) 


Reference 

Flow 

A1 

*2 

Adjustable  Parameters 

Phillips 
( 1977) 

Simple  Shear3 

(1  -  e)exp(-ets) 

f  exp(-gs) 

g  ,  a 

note  that 

r  =  -W  /V 

2.0X  1.0 

General 

(not  specified) 

(not  specified) 

Wagner 
( 1979) 

Simple  Shear3 

exp(-ns) 

0 

n,  a 

note  that  ^  ~  0 

General 

exp(-n  / 1  -  3) 

I  =  gl1  +  (1-<t)I2 

0 

Papana 

Stasiou 
et  al . ( 1983) 

Simple  Shear3 

a 

L  2 
(t  +  s 

0 

n,  o 

note  that  2  =  0 

General 

a 

0 

(n-3)  +  +  (1-«)I2 

aHere  s(t,t‘)  is  the  "magnitude  of  shear",  defined  for  v  =  y  (t)y,  v„  =  v_  =  0  as 
t,  X  yx  y  z 

s(t.t')  =  If  Y  (t" )dt" | . 

t  yX 

The  significance  of  the  factorizations  in  Eqs.  (1.8)  and  (1.9)  is  that  the  properties 
of,  say,  the  K-BKZ  model  now  depends  on  two  separate  functions:  the  linear  viscoelastic 
memory  function  f(t  -  t')  and  a  nonlinear  strain  function  W(I1(I2).  Standard  techniques 
exist  for  the  determination  of  M  from  measurements  performed  at  small  deformations  (Ferry 
1980).  In  this  report  we  will  demonstrate  how  Wd^,I2)  (or  respectively  A^d,j,I2)  and 
* Hl,I  ))  may  be  determined  without  the  restrictions  inherent  in  the  empiricism  in 
Table  1.2. 

Since  the  nonlinear  functions  to  be  determined,  W(I^,I2)  and  A^(Ij,I2),  i  =  1,2, 
are  functions  of  two  arguments  one  needs  experiments  that  cover  as  densely  as  possible  the 
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"obtainable  invariant  space"  given  by  all  possible  values  of  ( I  -j ,  1 2 )  subject  to  the 
incompressibility  constraint.  Such  experiments  may  presently  be  performed  in  just  one 
laboratory  namely  that  of  Meissner  at  the  ETH  in  Zurich  (see  Meissner,  Stephenson, 

Demarmels  and  Portman  (1982)).  The  experimental  method  can  be  thought  of  as  an  analog  for 
fluids  of  the  method  used  by  Rivlin  and  Saunders  (1951)  in  their  now  classical  experimental 
determination  of  the  strain  energy  function  for  natural  rubber. 


2.  KINEMATICS 


We  now  consider  the  class  of  flows  called  shear-free  flows.  We  use  the  notation 
proposed  by  Bird,  Armstrong  and  Hassager  (1985)  but  we  have  obtained  the  material  data  from 
Meissner,  Stephenson,  Demarmels  and  Portmann  (1982).  A  comparison  of  the  two  notations  is 
given  in  Table  2.1.  In  the  notation  of  Bird  et  al.  (1985),  all  shear-free  flows  may  be  put 
in  the  form: 

vx  =  -  —  £  (t )  ( 1  +b)x 

vy  =  “  2  "  b)V  (2.1) 

vz  =  £(t)z 

where  0  <  b  <  1  and  £(t)  is  any  function  of  time. 


Table  2.1:  Comparison  of  notation  in  shear-free  flows. 

b,  £:  parameters  of  Bird,  Armstrong  and  Hassager  (1985) 

m.  Eg!  parameters  of  Meissner,  Stephenson,  Demarmels  and  Portmann  (1982) 


Flow  Type 

b 

m 

Strain-rates 

Material 

Functions 

Uniaxial 

0 

1 

2 

•  • 

e  =  eo 

1  -+ 

1  ni 

biaxial 

0 

1 

£  =  -2iQ 

^1  = 

1  -+ 

3  ni 

plannar 

1 

0 

•  • 

6  “  "eo 

W1  = 

1  n’l 

W2  = 

K2 

We  have  analyzed  measurements  taken  by  Demarmels  (1983)  on  a  polyisobutylene  melt  in 


start-up  of  the  followinq  flows:  Uniaxial-,  biaxial-,  planar-  and  ellipsoidal-elongat ional 
flows.  Uniaxial  elongational  flow  is  obtained  by  choosing  b  *  0  and  £  positive,  and 
biaxial  elongational  flow  by  choosing  b  *  0 ,  but  £  negative.  Planar  elongational ,  also 

called  Mpure  shear  flow",  is  obtained  with  b  =  1  and  £  negative ,  and  ellipsoidal 

1 

elongational  flow  is  obtained  with  b  *  —  and  £  negative . 

We  wish  to  find  the  displacement  functions  for  start-up  of  the  shear-free  flow 

•  •  • 

described  in  Eg.  (1),  i.e.  we  imagine  that  £(t)  =  0  for  t  <  0  and  £(t)  =  £  a  constant 

for  t  >  0.  Then,  for  0  <  t'  <  t  we  get: 

,  -1/2  e ( 1+b) (t-t'  ) 
x  =  x  e 

,  -1/2  £ ( 1-b) (t-t  *  ) 

y  =  ye  2. 

,  +£(t-t') 
z  =  z  e 

and  for  *  <  t '  <  0  we  get 

,  -1/2  £ ( 1 +b ) t 
x  =  x  e 

,  -1/2  £ ( 1-b) t  ,, 

y  =  ye  (2.3) 

,  et 
z  =  z'e 

where  a  fluid  particle  P  with  coordinates  X£ ,  at  the  present  time  t,  had  coordinates 
xj  ,  i  =  1,2,3  at  a  past  time  t'.  We  define  s  =  t  -  t';  we  now  calculate  the 
displacement  gradient  tensor  for  0  <  t'  <  t.  In  the  notation  of  Bird,  Armstrong  and 
Hassoger  (19B5): 

f  e_1/2  e ( 1 +b ) s  0  0  1 


E(x,t,t* )  = 


-1/2  E ( 1-b) s 


and  we  find 


— £ ( 1 +b ) s 


E  •  F  =  B 


-E ( 1 -b )  s 
e  0 


'  V”  W'.'  V  V  V  r-  --  - 


We  now  find  the  first  invariant: 

X,  H  trig) 


Similarly  we  find: 


-e(l+h)s  -e ( 1-b)s  2es 
=  e  +  e  +  e 


1  /2  • 

e  '  e ( 1 +b ) s 


A( x, t , t 1  )  = 


1/2  £  ( 1  — b )  s 
e  0 


,+  .  -1 
A  «A  =  B 


£ ( 1+b) s  . 

e  0 


0  eE(1-b,S 


0  e 


The  notation  is  this  of  Bird,  Armstrong  and  Hassager  (1985), 
Similarly  we  find  the  second  invariant 


I.  =  tr ( B~ 1 ) 
2  “ 


E(1+b)s  ,  e(1-b)s  -2es 

=  e  +  e  +  e 


The  tensors  B  and  B  may  be  obtained  for  -»  <  t’  <0  merely  putting  "s"  equal  to 
"t"  in  Eqs.  2.5  and  2.8. 

For  the  incompressible  flows  defined  by  (2.1)  there  can  be  only  two  independent 
material  functions  in  start-up  flow: 


-+ 

n1  r 

-(  T 

-  T  )/t 

zz 

XX 

-  + 

• 

n  “ 

-(T 

-  T  )/E 

2 

ZZ 

yy 

and  in  fact  for  b  =  0  only  one  independent  material  function  can  be  found,  as 

-+  -+ 

n1  =  n2!  b  =  0  . 

This  will  become  obvious  when  we  write  down  the  expressions  for  n+ 


(2.11  ) 
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where  we  have  defined: 


*1  =  *1C11  +  ^2C1 2 


The  solution  to  Eq.  (16)  is  found  easily 


t  .  dn 

)  =  /  -  £  [ - Ids 

1  ^  G  ( s )  '■ds  ‘ 


(3.5) 


After  partial  integration  we  have 


•  -+ 

en. 


*  =  - - 

1  G(t) 


0  G(s) 


2  W  Vs 


(3.6) 


and  we  see  that  the  RHS  only  involves  the  experimental  data  (h^)  and  the  known  (measured) 
linear  viscoelastic  relaxations  modulus  G(t).  The  LHS  is  linear  in  the  strain 
functions.  If  we  use  the  same  procedure  for  the  second  material  function  we  find: 


•  -+ 

£n2 

*2  =  l*’lC21  +  *2C22  G( 


h  +  L 


^  £  r  dG  ( s )  >  -+ 

—  (“aHV8 


(3.7) 


0  G(s) 


where 


_  2£t  -£(1-b)t 


c21  =  e  -  e 


.  — 2et  e(l-b)t 


(3.8) 


c22  =  e  ”  e 


Eqs.  (3.6)  and  (3.7)  can  now  be  solved  explicitly  for  the  strain  functions  4>  i  and  $2' 
except  for  flows  with  b  =  0.  This  implies  that  it  is  only  possible  to  determine  a  linear 
combination  of  4^  and  4>2  if  we  only  have  data  from  uniaxial  or  biaxial  longational 
experiments . 

We  now  choose  a  specific  form  for  the  strain  functions;  which  is  the  same  form  as  used 
by  Bach  (1985)  to  obtain  accurate  approximations  to  the  Doi-Edwards  potential  function: 


>1  =  i  vw^i. 

m=1 


(3.9) 


*2  "  ^  VW^m 

m=1 


2- 


^TTTTT? 


vv 


7  •  ,  <  /i1;  t'v 


l  ; 

i 


r 


» 


VV’.-  V  V  VV\" 


'  V  J.-  !_'!• 


where  *  and  *2m  are  the  nodal  values  of  the  strain  functions,  and  Nm{I^,l2>  are  the 
global  shape  functions  for  the  isoparametric  4-node  element,  i.e.  the  strain  functions  are 
approximated  with  piecewise  bilinear  functions.  The  functions  are  defined  on  an  element 
mesh  covering  the  allowable  invariant  space  where  data  are  present.  We  may  now  write  eqs. 
(3.6)  and  (3.7) 


M  M 

K  vw^n  +  \  ww ci2  -  Hi(t) 


m=1 


m=1 


(3.10) 


M  M 

V  Nm^ 1 1 ,X2 )A1mc21  +  '  V X1 '*2 >*2mc22  *  H2(t) 

m=1  m=1 


where  Hi  and  H2  are  the  RHS  of  (3.6)  and  (3.7)  respectively.  The  arguments  of  the 

-e(1+b)t  ,  -f(l-b)t  .  +2et 


strain  functions  are  functions  of  elapsed  time  t  :  *  e 

e  ( 1 +b ) t  ^  e( 1-b)t  ^  -2et 


and  !•>  =  e 


+  e 


+  e  +  e 

+  e  .  we  can  calculate  the  Hi  and  H2  functions  at 


every  data  point  tn,  and  if  we  have  N  data  points  we  have  a  overdetermined  linear 
equation  system  with  2m  unknowns  and  2N  equations.  We  know  that  the  equations  (3.10) 
are  linearly  dependent  for  b  =  0,  and  one  of  these  may  a  priori  be  excluded  from  the 
final  equation  system.  The  least  squares  solution  for  the  2M  parameters  *im  can  now  be 
calculated  directly  using  any  standard  algorithm.  We  have  used  the  IMSL-routine  LLSQF.  In 
Appendix  we  give  a  short  description  of  the  FORTRAN  programs  we  have  developed  to  solve  the 
estimation  problem.  Here  we  only  note  that  it  is  convenient  to  use 


I*  =  log( I 1  -  2) 


and 


(3.11) 


I*  =  log(I2  -  2) 

as  the  independent  variables  for  the  strain  functions.  We  see  that  I*  and  1*  increase 
linearly  in  time  for  large  t,  which  gives  a  good  distribution  of  data  points  in  the 
(I*, I*)  space.  In  Figure  1  we  show  one  of  the  meshes  we  have  used  to  describe  the  strain 
functions.  Datapoints  are  indicated  by  (  +  ).  We  also  note  that  <t> ^  =  A^(I*,I*>  and 
*2  =  *2^r*'I2'  have  to  obey 

+l!0,0)  +  *2(0,0)  =  1  (3.12) 
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if  the  model  shall  conform  with  the  linear  viscoelastic  model  for  very  small 
deformations.  This  leaves  2N  -  1  parameters  to  be  estimated.  When  we  estimate  the 
parameters  in  eq.  (3.10),  with  the  restriction  (3.12),  we  find  a  very  high  (10^> 
condition  number  for  the  equation  system.  This  indicates  that  several  of  the  supplied 
equations  are  numerically  linearly  dependent.  In  the  two  flows,  uniaxial-  and  biaxial 
longational  (b  =  0),  we  know  that  it  is  impossible  to  determine  other  than  the 
combination: 

<(>  =  *,+  $2e’£t  (3.13) 

The  planar  data  do  apparently  not  supply  the  additional  information  needed.  If  the  mesh  is 
arranged  as  in  Figure  1  the  planar  data  do  not  supply  any  information  to  determine  the 
parameters  on  the  uniaxial  branch,  due  to  the  fact  that  the  shape  functions  belonging  to 
the  planar  branch  are  all  zero  at  the  uniaxial  branch.  If  we  avoid  this  situation  we  do 
find  a  decrease  in  the  condition  number,  but  only  a  factor  10,  we  may  conclude  that  it  is 
not  possible  to  obtain  two  numerically  reliable  and  Independent  strain  functions  with  the 
present  data  and  model.  To  obtain  numerically  reliable  parameters  we  simply  add 
constraints  to  the  strain  functions  by  requiring  equations  of  the  form 

-  af2  m  °r  lb  “  0)  (3.14) 

to  be  satisfied  in  least  square  sense  for  all  data  points  on  the  two  branches  corresponding 

to  b  -  0.  We  choose  a  to  be  5/2|  which  is  the  value  obtained  for  Poi-Edwards  strain 

functions  for  zero  deformation.  When  we  add  eq.  (3.14)  to  the  equation  system  we  obtain  a 

c  o 

very  large  decrease  in  the  condition  number,  from  10 J  to  10  .  This  Indicates  that  the 
calculated  parameters  are  numerically  reliable,  but  the  ultimate  test  of  a  model  is  of 
course  its  capability  to  reproduce  the  measured  material  functions.  The  obtained  strain 
functions  are  shown  on  Figure  3  to  Figure  6.  Doi-Edwards  strain  functions  are  shown  for 
comparison.  We  use  the  approximation  to  the  Doi-Edwards  strain  functions  obtained  by 
Currie  (1982).  Bach  (1985)  has  shown  that  the  approximation  do  not  give  deviations  in  the 
viscosity  functions  that  are  larger  than  5»,  and  this  accuracy  is  sufficient  for  our 
present  comparisons.  In  the  following  we  use  Curries  approximation  everywhere  the  Doi- 
Edwards  strain  functions  are  used.  The  most  striking  feature  of  the  calculated  strain 
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Fig.  3a 


optimal  strain  functions  in  comparison  to  the  Doi-Edwards  strain  functions  (+)  in  uniaxial 

•  “C  s 

;nsion  (b = 0,  c  >  0) :  a)  ♦  ;  b)  ♦  and  c)  $  =  ♦  +  $  *e  ;  where  the  s  corresponds 


Strain-dependent  memory  functions:  B  =0.000 
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Fig.  3c 
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Fig.  4a 


Figure  4.  The  optimal  strain  functions  (0)  in  comparison  to  the  Doi-Edwards  strain  functions  (+)  in  planar 
elongational  flow:  a)  4>  ,  b) 


Strain-dependent  memory  functions:  B  =1.000 
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Strain-dependent  memory  functions:  B  =0.333 


m 


& 

•H 

b. 


tO 

o 

to 

o 

to 

o 

to 

o 

CO 

CO 

CM 

CM 

«  < 

o 

o 

o 

o 

d 

o 

o 

o 

o 

d 

uoi-punj  uibj-^s 


-23- 


<PL0T79>  Release  1.6a  23-APR-84  08:41:43 


dependent  memory  functions 


Figure  6.  The  optimal  strain  function  (0)  in  comparison  to  the  Doi-Edwards  strain  function  (+)  in  biaxial 
elongational  flow  (b  =  0,  l  <  0):  a)  4>1>  b)  <j>2  and  c)  ^  +  ^e  £S,  where  the  s 
corresponds  to  the  acual  invariant  (see  eq.  (2.6)). 
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functions  is  that  they  are  not  "damping"-f unctions.  In  fact  they  increase  from  the  zero 
deformation  value  for  small  deformations.  This  is  true  for  both  strain  functions  in  all 
flows  considered.  It  is  also  obvious  that  the  Doi-Edwards  strain  function  would  give  a 
rather  bad  fit  to  the  data.  We  see,  for  example,  on  Figure  2  that  the  Doi-Edwards  strain 
function  is  typically  a  factor  2  smaller  than  the  optimal  strain  function  for  I*  larger 
than  2,  and  the  deviation  becomes  much  larger  as  I*  increases. 

The  ultimate  test  of  any  model  is  its  capability  to  reproduce  experimental  data.  We 
have  calculated  the  material  functions  in  all  flows  considered,  and  plotted  them  together 
with  the  data  on  the  Figures  7  to  10.  The  largest  derivation  is  found  where  the  data  shows 
the  largest  derivation  from  the  linear  viscoelastic  behavior,  which  is  in  uniaxial 
extension.  The  model  gives  a  reasonably  good  overall  fitt  especially  if  we  use  the 
"yardstick"  normally  used  when  model  predictions  and  data  for  polymer  results  are 
compared.  The  present  polymer  is  very  polydisperse,  with  a  weight  average  molecular 
weight  (f^)  that  is  2.14  times  the  number  average  molecular  weight  (Mn)  so  that 

-2.14  . 

We  may  also  notice  that  we  have  used  the  memory  function  obtained  by  Demarmels  (1983),  and 

his  procedure  for  calculation  of  this  function  may  not  be  the  best.  He  choses  the  time 

constants  a  priori,  and  then  he  calculates  the  coefficients  gk.  Whether  the  large 

deviation  found  by  Demarmels  (1983)  between  the  calculated  zero  shear  rate  function 

t/Tn 

h  +  (t)  “  'i  9kTk(l  “  e  )  and  the  measured  function  for  small  times  is  due  to  the 
estimation  procedure,  or  due  to  experimental  problems,  as  claimed  by  Demarmels  is  not 
clear.  It  is  not  likely  though  that  an  optimal  memory  function  would  imporve  the  models 
fit  to  the  present  data  significantly. 

We  may  conclude  that  it  is  not  possible  to  find  a  set  of  parameters  in  the  proposed 
model  so  it  Interpolates  the  present  data,  but  the  optimal  set  of  parameters  in  the  least 
square  sense  is  a  reasonable  approximation  to  data. 


0.0000  Gdo 


Comparison  of  model  predictions  (solid  curves)  and  data  (0  :  E  = 
=  0.02)  for  the  viscosity  function  in  uniaxial  elongational  flow 


Fig.  8a 


.018  and  :  e  =  -0.007)  in  planar  elonaational  flow:  a)  The  first  viscosity 


Fig.  9a 


Figure  9.  Comparison  of  model  predictions  (solid  curves)  and  data  (0  :  e  =  -0.122; 

X  :  e  =  -0.065  and  :  e  =  -0.012)  in  ellipsoidal  elongational  flow;  a)  The  first  viscosity 
function;  b)  The  second  viscosity  function. 
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Fig.  10 


re  10.  Comparison  of  model  predictions  (solid  curves)  for  the  viscosity  function,  and  data 
f  =  -0.104;  X  =  e  =  -0.043  and  :  e  =  -0.020)  in  biaxial  elongational  flow. 


4 .  ESTIMATION  OF  W. 

We  now  want  to  investiqate  the  consequences  of  an  assumption  of  the  existence  of  a 
potential  function  W  for  the  strain  functions,  i.e.  that 

A  -  A  =  -2-  (4.1  ) 

1  3t  5 

1  2 

The  W  function  has  to  be  at  least  C1  if  the  strain  functions  are  to  be  continuous. 
Alfeld  (1984)  has  developed  a  interpolant,  which  is  a  piecewise  cubic  function  defined 

on  a  mesh  of  triangles.  We  now  briefly  summarize  the  method  of  Alfeld  as  follows:  Each 
triangle  is  subdivided  into  three  subtriangles  (microtriangles),  see  Figure  11.  On  each 
microtriangle  a  cubic  function  is  defined. 


q(P)  = 


o+R+Y+^=3 


3!  „ 

q  !  R  !  Y  !  *  !  a«Y*b1b2b3b4 


P  =  )  b.V.  ; 

i=1  1  1 


is  the  point  of  evaluation.  The  "generalized  barycentric  coordinates"  bj^  satisfies 


i=i  '  ’ 


H  b.  =  0 
i=1  1 


The  "Bezier  ordinates"  cnpY*  are  functions  of  the  parameters  in  the  Clough-Tocher 
scheme.  The  Bezier  ordinates  are  given  explicitly  by  Alfeld  (1984).  The  parameters  are 
the  function  values  Q  and  the  directional  derivatives  at  all  vertex  nodes.  The 
directional  derivatives  j  ( cf )  are  defined  in  terms  of  the  partial  derivatives 

Dij(q)  -  '7q-eij  (4.5) 


where  e—  is  the  direction.  The  parameters  we  have  to  estimate  in  the  polymer  melt 


i 


it  is  one  third  of  the  macrotriangles 


4 


(1.1)  are  the  potential  function  values,  and  the  partial  derivatives  at  all  Conner 
vert  ices . 

To  calculate  the  strain  functions  everywhere  on  the  triangular  mesh  we  need  to 
calculate  the  partial  derivatives  of  (4.3)  with  respect  to  1^  and  I£.  We  find 


3q  _  y  3  !  3_  fvOnSwYiSi 

3x  .  ,  ,  ci !  6  !  Y  !  6  !  aSyS  3x  lbib2b3b4J 

a+8+y+6=3 


where  x  may  be  I*  or  I*  defined  in  Eg.  (3.11). 

It  is  trivial,  but  ouite  lengthy,  to  write  down  all  the  terms  that  evolves  from 
differentiation  of  the  products,  so  we  only  give  the  formula  for  the  partial  derivatives  of 


the  barycentric  coordinates  with  respect  to  the  invariants  I*,I£. 


The  coordina 


tes  b^  may  be  expressed  explicitly  in  I*, I*: 


b.  =  aQ.  +  au(I*  -  T* ( 4  ) )  +  a2.(I*  -  I*(4))f  i  =  1,2,3 


where  ag^  =  —  for  all  i,  and 


*11  -  k  (I2(2)  -  I2<3));  a2 1  =  k  (I1(3)  -  I*(2)> 

a12  -  k  (I2(3)  -  I2(1>)'-  a22  =  k  (I1(1)  -  I*(3,) 

a  1 3  =  k  (I2(1)  "  I2(2));  a23  =  k  (I1(2)  *  I1(1)) 


See  Fiqure  11  for  a  geometrical  interpretation  of  above  formula 

(I*(i),I*(i)  )  =  =  vertex  coordinates  of  microtriangle 

We  now  have  the  simple  result: 


3!*  =  aU  an<1  31*  “  a2i 


We  may  now  formulate  the  estimation  problem: 
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where  the  3M  parameters  are: 


=  r rw( P  ,r*iU  );...,fW(P  > n 

i  '  a i*  p  '  ax*  p  '**•'  1  m  si»  p  '  ai*  p 

112  1  1  M  2  M 


=  f  ( °’  -j , '  '"”2 !  '"'3  ^ 


where  M  is  the  number  of  vertices.  We  note  that  W  does  not  enter  into  the  objective 
function  (4.8),  and  consequently  we  must  supply  at  least  one  w  value.  This  may  be  picked 
arbitrarily.  We  take 

W(I*,I*)  =  0  at  (I*, I*)  =  (0,0)  . 

We  want  to  approximate  the  potential  function  W  by  (4.3): 


W(I*,I*)  = 


«+4+y+*=3 


3 !  ™  R  y  * 

a!«!Yl*  !  Cn#Y*  b1b2b3b4 


(4.10) 


where 

Coi8Y*  =  ^Y*1"1 

(i*,1*)  =  b,vt  +  b2v2  +  b3v3  +  b4v4 

and  are  the  vertex  coordinates,  including  the  centride  (V4),  and  b^  are  the 

Barycentric  coordinates  of  the  triangle  containing 

If  the  model  shall  conform  with  the  linear  viscoelastic  model  for  very  small 
deformations  we  must  require: 


f2«_i  +  f2iU 

5I  aj 

1  2 


<11'I2)  = 


This  leaves  3M  -  2  parameters  left  to  be  estimated.  The  objective  function  is  a  sum  of 
squares  of  cubic  functions  and  the  optimization  problem  is  consequently  supposed  to  be 
simple,  and  we  may  choose  any  standard  unconstrained  algorithm,  in  particular  one  that  does 
not  require  analytical  derivatives  of  the  objective  function  with  respect  to  w.  These 
derivatives  are  rather  lengthy  to  write  down. 
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In  Appendix  1  we  give  a  short  description  of  the  program  developed  to  estimate  the 


parameters  m  in  (4.9). 

On  the  Figures  13  to  16  we  present  the  calculated  potential  function,  and  its  partial 


and 

z2“ 

We 

have  used 

3  W 

3w 

1  3w 

3h  = 

3I1 

h 

-  2  ”  31* 

3w 

3W 

1  3w 

3I2  = 

3I2 

1  2 

-  2  _  31* 

-I! 


(4.11) 


-T* 

"2 


(4.12) 


The  invariant  space  has  been  transformed  into  a  regular  rectangular  mesh,  as  required  by 
the  plotting  routines.  We  use  a  dimensionless  time  T  =  2et  as  the  first  independent 
parameter,  and  a  "flow  type"  parameter  £  as  the  other: 

+1  biaxial  longational  flow 
0  planar  longational  flow 
-1  uniaxial  longational  flow 
In  Figure  12  we  have  made  a  draft  of  the  transformation,  and  the  two  different  view- 
directions  used.  The  zero  deformation  point  is  transformed  into  a  line.  The  maximum  t 
value  is  approximately  6,  which  corresponds  to  the  largest  invariants  in  the  data  (see 


Figure  1).  We  also  note  that  the  different  figures  have  different  scales  on  the  "z-axis", 
r  3  w  i  ^  ,  ,  f  3w 


in  fact  =  and  (gy- )g  =  0.34. 

We  initially  tested  the  influence  of  the  number  of  polynomials  used  to  describe  W. 


We  found  no  improvement  using  more  than  3  cubic  polynomials.  (We  tried  to  use  up  to  33 
cubic  polynomials). 


3W 


The  most  characteristic  feature  of  the  jrz - function  is  that  in  uniaxial 

3Ii 


elongafional  flow  it  increases  rapidly  above  the  zero  deformation  value  for  moderately 


all  t,  and  for  large  t  it  decreases.  The  opposite  is  true  for  biaxial  elongational 


3w 


flow,  and  here  we  see  a  monotanic  decrease  in  time.  The  - - function  increases 

3 1 2 


initially  for  all  types  of  flow. 


Fig.  13 


Figure  13.  The  potential  function  W  is  transformed  into  the  (t,F)  plane  and  plotted  as 
surface.  The  view  direction  is  Northwest  (NW)  (see  Figure  12  for  definition). 


The  derivatives  are  plotted  in  the  same  way  as  w  on  Figure  12: 


Fig.  15 


The  W-function  is  plotted  with  view  direction  Southeast  (SE) 


re  16.  Tt>e  derivatives  are  plotted  as  W  in  Ficiure  15: 


To  obtain  an  idea  about  how  the  strain  function  depends  on  the  invariants  we  have 
3  w  3w 

plotted  the  yT‘  t  and  the  ir - function  on  Figures  17  to  19  as  functions  of  Ii  with 

dI1  3I2 

b  as  discrete  parameter.  We  have  also  plotted  the  "effective"  strain  function  in  the  two 
extreme  flows  with  b  =  0.  The  Doi-Edwards  strain  functions  have  been  plotted  for 
comparison.  We  notice  again  the  initial  increase  in  ( )  in  uniaxial  flow.  In  Figure 
17c  the  effective  strain  function  combination  is  plotted,  and  we  notice  a  good  coincidence 
between  the  estimated  strain  function  and  Doi-Edwards  for  very  small  deformations  but  very 

soon  the  derivations  become  very  large  (>100%).  There  is  no  overshoot  in  planar  extension 

„  3w  ,  „  3W 

for  tt— #  but  now  there  is  for  m: — • 

31,  il 

3  W 

There  is  a  very  small  overshoot  in  -r—  in  biaxial  extension.  (It  is  too  small  to  be 

0  I 

seen  on  the  3-D  plots).  The  effective  strain  function  shows  a  larger  overshoot,  as  it  is 
3  W 

now  the  - -  function  that  is  the  important  strain  function.  Again  there  is  a  large 

3I2 

deviation  between  the  estimated  strain  function  and  the  Doi-Edwards  function. 

Figure  20  to  23  shows  how  the  estimated  potential  function  describes  data,  both  in  a 
linear  scale  and  in  a  logarithmic  scale.  The  latter  is  often  used  as  the  customarily  large 
deviations  between  data  and  model  predictions  look  quite  small  on  such  plots. 

The  deviations  in  uniaxial  flow  are  approximately  the  same  as  we  observed  for  the  more 
general  model.  In  planar  extension  we  see  that  the  curve  for  the  largest  E  -  value  has 
improved,  but  the  curves  for  the  two  smallest  E  -values  have  become  worse.  The  data, 
especially  for  the  smallest  e  -  value,  do  not  seem  to  be  accurate.  The  second  viscosity 
function  has  improved  for  all  E  -  values.  We  may  notice  that,  although  there  is  a  large 


negative  deviation  from  the  linear  viscoelastic  behaviour  for  the  largest  e  -  value,  the 
model  is  apparently  much  better  to  describe  this  kind  of  deviation.  Again  in  ellipsoidal 
extension  we  see  a  small  improvement  for  the  largest  E  -  value,  and  an  overall  good  fit 
for  both  viscosity  functions.  In  biaxial  extension  we  observe  the  same  deviations  as  we 
saw  when  we  did  not  assume  the  existence  of  a  potential  function. 

We  may  now  conclude  that  a  potential  function  "exists",  in  the  sense  that  we  obtain  a 
fit  to  data  which  is  as  good  as  if  we  do  not  assume  the  existence.  In  fact  we  get  some 


-<16- 


Strain  function  B  =  0.000 


Figure  17.  The  strain  functions  are  plotted  as  functions  of  I1  in  uniaxial  flow: 


Fig.  17b 
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Fig.  I7c 
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Strain  function  B  =  0.000 


The  strain  functions  are  plotted  as  functions  of  1 1  in  biaxial  flow: 
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function  in  linear  scale,  b)  viscosity  function  in  logarithmic 


Fig.  20b 
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.065  and  :  F  =  -0.012)  in  ellipsoidal  elongational  flow:  a),  c)  first  viscosity 
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re  23.  Comparison  of  model  predictions  (model  (1.1))  and  data  (0:6=  -0.104, 

c  =  -0.043  and  :  c  =  -0.020)  in  biaxial  elongational  flow.  a)  first  viscosity  function  on 
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function . 


5.  DISCUSSION 

Our  investigations  of  the  factorized  single  integral  constitutive  equation  suggests 
that  the  assumption  of  existence  of  a  potential  function  for  the  strain  functions  does  not 
limit  the  model  in  its  ability  to  describe  experiments.  We  find  equally  good  descriptions 
of  the  general  elongational  data  for  the  investigated  PIB-polymer  melt,  with  the  assumption 
of  a  potential  (Eq.  1.1)  as  without  (Eq.  1.2). 

We  also  found  that  with  the  factorized  single  integral  constitutive  equations  and  the 
linear  viscoelastic  memory  function  given  by  Demarmels  (1983)  it  is  iiqpossible  to  adjust 
the  strain  functions  to  describe  the  experiments  of  Demarmels  (1983)  without  significant 
systematic  deviations.  There  may  be  several  reasons  for  this: 

i)  The  linear  viscoelastic  memory  function  given  by  Demarmels  may  not  be  the  optimal 
function.  For  example  we  do  not  feel  comfortable  with  the  oscillatory 
"measurements”  in  the  frequency  range  from  10  s  and  10  s  artificially 
constructed  from  time  temperature  superposition.  Also  the  shortest  relaxation 
time,  10"3  sec.  should  be  compared  with  the  fact  that  a  constant  strain  rate  is 
not  established  in  the  experiments  before  after  about  3  sec.  It  might  have  been 
preferable  to  represent  relaxation  times  so  short  compared  to  typical  times  in 
the  experiment  by  just  a  Newtonian  contribution  and  estimate  this  contribution 
directly  from  the  elongational  flow  experiment. 

ii)  One  can  not  rule  out  systematic  errors  in  the  experimental  data.  In  particular 
what  we  fit  is  not  raw  data,  but  rather  the  results  of  some  analysis  applied  to 
the  measurements.  This  analysis  is  often  different  in  differnt  flow  situations 
even  for  the  same  apparatus  (Demarmels  1983,  Ch.  5)  and  consequently  a 
possibility  for  different  types  of  systematic  errors  in  the  material  functions 
arises. 

iii)  The  most  interesting  possibility,  and  the  one  we  consider  most  likely  is  that  the 
factorized  single  integral  models  with  a  memory  function  determined  entirely  from 
linear  viscoelastic  measurements,  simply  does  not  have  sufficient  flexibility  to 
represent  nonlinear  material  behavior  over  the  entire  invariant  space.  We 
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believe  it  is  the  first  time  that  possibility  has  been  so  clearly  established  on 
the  basis  of  experiments  with  "ever  increasing  deformations.”  The  factorized 
single  integral  models  are  known  (Demarmels  1983)  to  be  unable  to  describe  flows 
with  sudden  changes  in  flow  type. 

The  strain  functions  calculated  from  the  estimated  potential  function  differ  greatly 
from  previously  proposed  strain  functions;  see  Tables  1.1  and  1.2.  The  optimal  strain 
functions  for  the  particular  PIB  melt  have  maxima,  with  values  well  above  1,  so  the  concept 
of  the  strain  functions  as  being  "damping-function"  is  not  true  for  the  particular  model, 
and  material.  This  point  was  noted  also  by  Demarmels  (1983).  The  general  elongational 
flow  data  used  here  have  the  unique  feature  that  they  cover  major  parts  of  the  invariant 
space.  In  contrast,  the  more  common  transient  shear  experiments  are  limited  to  1^  =  I2. 
The  general  elongational  experiments  have  two  disadvantages,  however.  First  it  is  not 
possible  to  obtain  data  for  times  that  are  even  comparable  to  the  smallest  time  constants 
of  the  fluid.  The  present  rheometer  can  not  give  measurements  with  constant  £  for  times 
smaller  than  approximately  3  sec.,  and  the  smallest  time  constant  is  0.001  sec.  Second  it 
is  not  possible  to  obtain  measurements  for  large  e  -  values,  or  invariants. 

If  one  were  to  search  for  a  more  general  model  our  data  comparison  indicates  that  one 
should  build  on  Eq.  (1.1).  We  point  out  two  possible  generalizations: 

i) 

t  aw.  aw. 

r  ( t )  =  -n  Y  +  f  )  M ,  (t  -  t 1  y...  +  -=rr-  Y  ldf  (5.1) 

~  s~  .1  ai  ~[0]  ai  ~ 

1  1  2 


where 


M^s) 


n . 

— j  exp( -s/X . ) 

X f  1 

1 


(5.2) 


and  v  is  the  rate-of-strain  tensor.  Here  the  short  relaxation  times  have  been 
"strangled"  and  represented  by  a  Newtonian  term.  In  other  words  to  a  first  approximation 
^  where  the  summation  is  on  the  relaxation  times  small  compared  to  the  shortest 


time  in  the  experiment,  here  about  3  sec 


ii)  One  could  add  a  term  similar  to  the  Curtiss-Bird  (1981)  e-term.  This  involves 
the  rate-of-strain  tensor  contracted  twice  with  a  fourth  order  tensor,  and  it  is  not 
immediately  obvious  how  one  would  determine  this  fourth  order  tensor.  The  method  proposed 
by  Currie  (1982)  involves  second  derivatives  of  a  potential  function.  Alfeld  (1984b)  has 
derived  an  explicit  C2-Clough-Tocher  scheme  so  this  approach  is  feasible,  but  a  serious 
question  is  if  it  is  possible  to  implement  models  with  this  complexity  in  general  flow 
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APPENDIX 


A  brief  documentation 

FORTRAN  IV  programs  for  the  use  in  parameter  estimation  in  a  Kaye-BKZ  polymer 
melt  model.  The  theory  and  sample  results  for  a  specific  PIB-polyrner  is  given  by  Bach 
and  Hassager  (1984a). 


Estimation  programs. 

The  programs  ESTIM  ,  and  ESTW  are  the  main  programs  for  estimation  of  two 
independent  strainfunctions,  and  for  estimation  of  a  potential  function  respectively.  The 
input  is  raw  start-up  data  from  any  elongational  flow.  The  time  dependent  behaviour  is 
in  the  model  assumed  to  be  described  by  the  linear  viscoelastic  memoryfunction,  which  is 
supplied  as  a  FUNCTION.  Only  a  very  little  programming  effort  is  required  to  make  the 
estimation  programs  take  other  material  data.  e.g.  shear  data  (stress  relaxation,  stress 
growth  experiments). 

ESTIM  estimates  two  independent  strainfunctions  4>j  and  These  functions  are 
defined  by  the  shape  functions  belonging  to  an  element  mesh  covering  the  invariant  space. 
ESTIM  uses  the  functions  GPIB  and  FPIB.  These  are  the  relaxations  modulus  and  the 
memory  function  respectively.  When  the  linear  least  squares  equation  system  is  set  up  by 
ESTIM  the  IMSL  routine  LLSQF  is  called  and  finally  the  solution  is  written  onto  a  file. 

A  number  of  routines  from  the  Newtonian  Flow  program  package  NEWTON  (Bach 
and  Hassager  (1984b))  have  been  used  in  ESTIM. 

ESTW  estimates  the  parameters  in  a  C1  Clough-Tocher  scheme  Alfeld  (1984)  for 
the  potential  function  W.  The  Clough-Tocher  scheme  is  implemented  in  the  subroutine 
ClIN’T.  This  routine  calculates  the  partial  derivatives  and  jjj;  ,  and  optional  the 
potential  function  value,  at  any  point  in  the  covered  invariant  space.  The  invariant  space 
is  transformed  to  obtain  a  better  distribution  of  data  points.  We  use  7"  =  log(I  -  2)  and 
II'  =  log(II  -  2).  The  partial  derivatives  are  related  by  : 


cBV  dW  1 _  dW 

di  "  di-  (7  -2)  “  dfe 

aw  _  aw  i  aw 
dll  ~  dll ■  (//  2)  =  dTTe 

It  would  in  principle  be  possible  to  use  the  same  method  as  in  ESTIM  but  as  the 
equations  system  would  be  very  difficult  to  set  up.  a  straight  forward  minimazation  method 
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is  used  instead.  The  1MSL  routine  ZXSSQ  performs  very  well  on  the  actual  problem.  The 
residuals  are  calculated  by  the  FUNCTION  FC2. 

All  subroutines  used  are  described  in  the  following.  They  are  listed  in  alphabetic 
order.  All  parameters  in  boldface  are  to  be  set  by  the  user,  and  all  underlined  parameters 
are  to  be  reset  by  the  user,  as  they  are  changed  by  the  routine. 

ClINT(X,  Y,  IOPR,  MP.  W,DW,  DDWI.DDW2,  DDW 3,  F,  DFX,  DFY,  1ER) 

ClI.NT  calculates  the  function  value  F  (only  for  IOFR  —  0  ),  and  the  partial  derivatives 

-  DFX,  and  §£  =•  DFY  of  the  interpolant  at  (X,Y).  If  IOPR  >  0  then  ClINT  assumes 
that  (X.Y)  is  in  element  number  —  IOPR.  If  IOPR  <  0  then  ClINT  searches  all  elements 
until  the  one  containing  (X.Y')  is  found.  If  no  element  is  found  IER  equal  to  1  will  be 
returned,  else  IER  —  0  is  returned,  and  IOPR  is  set  to  the  actual  element  number  found. 
Note  that  this  option  can  be  used  to  enhance  the  effectiveness  greatly  in  situations  where 
ClINT  is  to  be  called  with  the  same  set  of  (X,Y)  many  times.  Topology  data  is  supplied 
in  the  Cl  DATA  -  COMMON  block,  and  the  data  to  be  used  in  the  interpolation  scheme 
are  suplied  in  W(MP),  and  DW(MP,2)  .  The  arrays  DDWl(MP,2),  DDW2(MP,2),  and 
DD\V3(MP.2)  are  the  directional  derivatives  at  the  local  nodes  1,2  and  3  respectively. 


FC2(OMEG  A,  NRES,  KPAR, RESV) 

FC2  is  the  EXTERNAL  routine  required  by  ZXSSQ.  The  parameters  are  suplied 
in  OMEGA(NPAR),  and  all  residuals  are  to  be  returned  in  RESV(NRES)  for  any  given 
parameter  vector  suplied  by  ZXSSQ.  Every  datapoint  contribute  with  two  residuals,  except 
for  flows  with  B  0. 

FMEMOfS.  Td,  GO) 

FUNCTION  FMEMO  calculates  the  Doi- Edwards  memory  function.  S  is  the  inde¬ 
pendent  argument  of  the  memory  function,  and  Td  :  -  Td  is  the  disengagement  time  .  and 
CO  :  Go  is  the  elastic  modulus.  The  function  value  is  returned  in  FMEMO. 

FPIB(S.rrf.GO) 

FUNCTION  FPIB  calculates  the  memory  function  obtained  by  Demarmels  (1983).  S 
is  the  argument  of  the  memory  function,  and  here  is  GO  the  zero  shear  rate  viscocity  ijo- 
Td  is  not  used.  The  function  value  is  returned  in  FPIB. 

G AUSSW(N, ROOT.  W  EIGHT) 

SUBROUTINE  GAUSSW  calculates  the  N’th  Gauss-Legendre  rule;  that  is  the  abcis- 
sae:  ROOT(N)  in  the  interval  1  ,  and  the  corresponding  weights:  WEIGHT(N). 


GMEMO(S,  TD,  GO) 

FUNCTION  GMEMO  calculates  the  Doi-Edwards  relaxation  modulus. (See  FMEMO 
for  the  explanation  of  the  arguments). 

GPIQ(S,Td,GO) 

FUNCTION  GPIB  calculates  the  relaxation  modulus  obtained  by  Demarmels  (1983) 
for  a  commercial  PIB  polymer.  (See  FPIB  for  the  explanation  of  the  arguments). 

LAMDA(I,II,  L) 

SUBROUTINE  LAMDA  calculates  the  eigenvalues  in  any  flow  using  the  analytical 
formulae  derived  by  Currie  (1982).  I,  II  are  the  first  and  second  invariants  respectively. 
Note  Curries  definitions  !  It  shall  be  noted  that  Curries  formulae  are  not  suited  for  nu¬ 
merical  evaluation  for  large  invariants.  Numerically  it  would  be  better  to  solve  the  cubic 
equation  for  the  eigenvalues  directly  Bach  (1985). 

UPOT(l,II,N,I,F/l,E/2) 

FUNCTION  UPOT  evaluates  the  Doi-Edwards  potential  function  W  using  a  N’th 
degree  Gauss-Legendre  quadrature  formula.  The  partial  derivatives  of  W  with  respect  to 
the  first  (I)  and  second  (II)  invariant  are  calculated  using  the  analytical  formulae  derived 
by  Bach  (1985).  The  eigenvalues  L(3)  are  calculated  by  LAMDA.  The  partial  derivatives 
~W1  an<^  §71  at  (MI)  are  returned  in  Fll  and  FI2  respectively.  The  value  of  the  W  - 
function  at  (I, II)  is  returned  in  UPOT. 


Utility  programs. 

The  programs  START  and  STARTl  are  the  interactive  main  programs  used  to  calcu¬ 
late  material  functions  and  to  test  the  obtained  models. 

START:  Interactive  MAIN  program  . 

START  calculates  start-up  material  functions  in  shear  flows,  and  in  all  kinds  of  elon- 
gational  flows  for  1)  The  exact  Doi-Edwards  model,  2)  Curries  (1982)  approximation  to 
Doi-Edwards  model.  3)  General  Kaye-BKZ  models,  with  two  independent  strainfunc- 
tions  given  as  nodalpoint  values  on  an  elementmesh.  The  exact  Doi-Edwards  potential 
function  is  calculated  by  UPOT.  The  memory  integral  is  evaluated  by  a  very  cautious 
Romberg  routine  called  DCADRE  (IMSL).  START  writes  the  necessary  information  to 
the  PLOT79-interface  routines  (GPLT,  GPLTL,  GPLTLL).  START  can  also  be  used  to 
test  model  prediction  against  viscometric  data. 
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START1  :  Interactive  MAIN  program. 

STARTl  calculates  the  start-up  material  functions  in  shear  and  all  kinds  of  elonga- 
tional  flows  for  the  Kaye-BKZ  polymer  melt  model  with  a  potential  function  W.  W  ,  and 
its  derivatives,  are  calculated  be  C1INT.  The  memory  integral  is  evaluated  by  DCADRE 
( IMSL-routine).  STARTl  writes  the  necessary  information  to  the  PLOT79-interface  rou¬ 
tines  (see  START).  STARTl  can  also  be  used  to  test  model  predictions  against  viscometric 
data. 
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